{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:44.301431Z",
     "iopub.status.busy": "2022-01-02T18:59:44.300782Z",
     "iopub.status.idle": "2022-01-02T18:59:44.303686Z",
     "shell.execute_reply": "2022-01-02T18:59:44.303074Z",
     "shell.execute_reply.started": "2021-12-10T11:52:23.421825Z"
    },
    "papermill": {
     "duration": 0.022357,
     "end_time": "2022-01-02T18:59:44.303801",
     "exception": false,
     "start_time": "2022-01-02T18:59:44.281444",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.011631,
     "end_time": "2022-01-02T18:59:44.327795",
     "exception": false,
     "start_time": "2022-01-02T18:59:44.316164",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## The original data\n",
    "\n",
    "We first read original data\n",
    "\n",
    "(The data used in this notebook have been uploaded here https://www.kaggle.com/elmahy/pems-dataset)\n",
    "\n",
    "Its shape is # (sequence_length, num_of_vertices, num_of_features) (16992, 307, 3)\n",
    "\n",
    "We have 307 detectors each detector detect three features such as speed of cars, congestion,etc. \n",
    "And we store that data every five minutes. \n",
    "In 5 mins, our data will be (1, 307, 3)\n",
    "In 1 hour ,  our data will be (12, 307, 3)\n",
    "In 1 day ,our data will be (12 * 24, 307, 3)\n",
    "In 59 days, our data will be (12 * 24 * 59, 307, 3) =  (16992, 307, 3)\n"
   ]
  },
  {
   "attachments": {
    "66f43f1d-1b32-4eb8-8c59-814f64b9d4ae.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAAB/CAIAAACmDsuiAAAAA3NCSVQICAjb4U/gAAAAGXRFWHRTb2Z0d2FyZQBnbm9tZS1zY3JlZW5zaG907wO/PgAAIABJREFUeJzsnWdcVEcXh8/2DiydpfcqxY69xa4RE4NRY29R35ioMYnGSKyxJWoSoyRGY+xiSVCwoRgBFakL0ntZ2i6wsJUt9/1wcYOAikoxMM8HfnfnzpyZO7vc//RDwDAMEAgEAoHoSRC7ugAIBAKBQHQ2SPwQCAQC0eNA4odAIBCIHgcSPwQCgUD0OJD4IRAIBKLHgcQPgUAgED0OJH4IBAKB6HEg8UMgEAhEjwOJHwKBQCB6HEj8EAgEAtHjQOKHQCAQiB4HEj8EAoFA9DiQ+CEQCASix4HED4FAIBA9DiR+CAQCgehxIPFDIBAIRI8DiR8CgUAg3lJmzpxJIBDmzJnT7paR+CEQCES3Ra1W79y508PDg0aj6evr9+vXj8/nA0CfPn0WLFgQGBjo5OREJBIJBMKXX36pS1VVVWVlZUUgEAgEQk5ODgCEhoYSCIR9+/aNHDnS3NycRqNZW1svWLCgtLRUKBQuWrTI1NSURqN5eHj89ttvmZmZEydONDc3xy1cv34dACIjI0ePHm1sbEyj0Xg83syZMwsLCwFg2rRphOdQUlLScTVD7jjTCAQCgehCMAybMWPGlStXuFzunDlzmEwmn8+vqqoqLCxMSEj45ptvli9f7ufnJxaLhUKhLpVWq50zZ05FRUVTU5cuXerTp8/9+/eTk5P9/f3JZHJYWNjx48fj4+PJZHJiYqKtre2QIUNCQ0OXLFnyv//9Lzs729/f/8qVK3jympqaCRMmKBSK6dOne3t7nzhx4ty5c9XV1Tdv3nR1dR0wYAAAlJSUlJaWEonEfv364aloNNqbPLtWqyWRSC+KgUAgEIjux7Vr1wCAw+EUFhY2Df/+++85HI5CocA/9unTBwC++OIL/GNQUBCBQNizZw+uEdnZ2SqVysjIaPv27ffv35dKpXi03bt3N5WSvLw8DMPw7qOpqalGo6mvr8dvhYeHJycn49e///57amrqlClTAGDWrFlNS7Vx40YAYLFYTQMDAwMBIDAwcM2aNVwu18TEZP369Vu3bgUALy8vPM7q1asBYPDgwXgWv/zyi7u7O4lESktLe0HloGFPBAKB6J6Eh4cDgLGxcWBgIIvFcnBw2Lp1q0ajuXTp0qRJk1rtV0VERGzZsmXTpk3jx4/XBd67d08kEk2fPn3IkCFMJhMPlEqlAECn0wHA3Nzc3t4eAAYNGgQAlZWV+GCpjl69egUEBADAwoULvby8QkND33///d9//72NDxISEhIWFtarV6+qqqrdu3dnZWUBAIFAaDXyqlWrHBwcpkyZQqFQXmATiR8CgUB0T8rLywEgPz+fQqEsWLCgoqLim2+++fLLL2NiYqZPn94yvkAgmDVr1ujRozdv3tw0/NKlSx4eHm5ubrqQkJCQHTt2UCiUoUOHAgCbzcbDORwOflFTU9PUAoFAmDt3Lo/Hmz9//u+//+7i4hISEvLdd9+18UGsrKySkpLu3bvn6uoKAKWlpQBAJLauX999993Vq1cvX77s5OT0AptI/BAIBKJ7gmsSkUgMCwv76aefli5dCgCnTp2iUqkTJ05sGT80NLSyslIikUydOnXlypV44JIlS86cOaMTS61Wu2nTpg8++IDFYoWHh3t7ewOARCLB7+qGOrlcblPLDx8+DAgIkMlkR48eXbBgQVBQEADs3r0bw7C2PEj//v3xfqqZmRkAKJXKpnfVanXTj+PGjWuLTbTgBYHo5jx48GDTpk1ubm5kMjklJaV///719fXp6ek7duzAFxoguisDBw48fvw4AODrPvC/Mpls3LhxLBarZXxcih48eNA0MDIyEgDee+89ABCLxbNnz7527Zqnp+eVK1ecnJzEYvG+ffvKy8sLCgrs7OzwtKampk5OTjKZTGcEn/NTqVQqlYpGo+G3ntd1awmVSsUv8KFOfDxT17lMTExsGpnBYOAXdXV1f//9NwBMmzZN1zfVgcQPgejmJCUlTZkyZfXq1aGhoampqTt37sQwbOTIkfgIEqIbM3v27F27duXn5wcEBPj6+gYHBwOARCLRdeOWLFmiUqkKCgoA4Nq1a3369Pnhhx8+/fRTAEhNTe3VqxcALF68OCIiwtfXFwAWLVqEL6Jhs9nLly8HAK1W6+bmlpGRMXr06L59+16+fBkAli1btnDhQl2HbN++fSwWi0wmS6XSSZMmDRw48LfffgOAqVOnPm/e7sWYm5sDQFFR0YIFC+rq6mJiYlqNJhAIPvroIwDIzs5uOQSKxA+B6ObY2dkNHz4cAPh8vpeXFwAQCIRZs2YZGBh0ddEQHQubzb537966detu3boVFRXl4uLSv3//o0eP4ostAeDPP//UDSGmpqampqYKhUJc/HTcvHlzxowZ+LVuePPRo0e6CGFhYefPn8en2RwdHT/99FMfHx9/f39dhNu3bwPA4sWLc3Nzk5OTo6KiLC0t165d22xmse1YWlpu3br1xx9/vHr16pgxY2bPnn3q1KlXNUJo45Ar4k1YsmRJRUWFu7t7Wloak8m0t7d/8uSJmZkZ3vxBIDqHmTNnjhkzZvHixV1dEESXERAQIJVKb9682cb4ycnJvr6+0dHR+DLO7gTq+XU4SqWyuLj45MmTxsbGkydPDggImDlz5tmzZ/GhBgSic5DJZHl5efhAFqLH4u/vP3DgwLbHVyqVW7dubdqN6zYg8etwJBLJxo0bjY2NhUJheXk5Pu7k7u7+4mW4CET7kpaWptVqPT09W97KzMxMS0vr16+flZVV5xcM0ZmsX7/+leL379+/f//+HVSYrgVtdehwjIyM8K0wKSkpFAoFX2Xg4+PTt2/fri4aogfB5/Pt7e1brnkDAEdHx0OHDjU0NHR+qRCIrgL1/DoPPp/v7u7+4kMHEIgOQrd4ryl5eXn//POPkZGRRqPBT+hAIHoIqOfXeaSkpKAZF0SXoNVq+Xx+0xM6AOCff/45cODA1KlTc3NzPT09X2/ROQLxHwWJXyehVqszMjKQ+CG6hB9//LG2tjY2Nla371itVn/33XdfffWVoaEhnU738fHp2hIiOhOFQvH55587OjoyGAx9fX0/P7/Dhw83ixMSEtK7d286nW5oaBgYGFhcXNzSdVGzJEeOHPHz89PT02MymU5OTmvXrlUoFLq7JSUlLZ0WLV68uFkIfvhLJ4CGPTuJnJwchUKBr3YBgPT0dKFQqJsLtLOzYzAYUVFRVlZWaCEMot1ZvXo1fvK9jszMTCaTiW8Wfvjw4ccff9xFRUN0AZs3b967d6+Zmdm6detKSkqOHz/+8ccfe3t76/YzhISE4Hv7Jk+enJOTc/78+djYWH19/eTk5Kaui5RKpe4UtPDw8OXLl5NIpJUrV9Lp9J9//vn77783MjLasGFD06yZTGbLtaNeXl74uWUA4ODg0LEPr+P1PGUgXpXTp0+PGjVK97GkpKRPnz4qlQrDsHfffffBgwdarXb9+vVnzpzpujIiehCxsbG4Q5m7d+8OHTpULBYLBIKuLhSik8AdLAwcOJDP5+MntpDJZPyUzo0bN2IYho9RzZs3D8Owmpoa3HUDTkvXRbjNAwcOAIC+vn5UVFRcXJytrS0AfPLJJ3gqlUpVXFwMAK6urk1LsmjRIgBo9b0nlUqHDx9ubm5OpVLpdLqvr++JEydwD0czZsxYsmQJi8WytbUNCQm5fv26q6sri8WaOHFiRUVFGysBDXt2Bjk5OefOnZPL5Q8fPsRDLC0tORxOdXV1amqqWq2uqqoCAJFINHXq1C4tKaKn4OnpqdFo1q9fT6VSfXx8du7cKZfLu7pQiE5i9erVenp6Dx8+9Pb2njRpEpfLDQ8P1x2JKZFIUlJSAADvohkYGHh4eOC3XuC6KDAw0NXVVSwWDxkypG/fvsXFxbt3754wYUKzrPPy8jgcDpfLHTJkCK67ALB8+XI6ne7g4LB69erq6mo8sKGhIScnJyAg4Jtvvnn33XeTkpLmz5+P37148WJ8fLyzs3NhYeGcOXNmzJjh4uKC+9fdtGlTW2vhDVoPiDfiww8/TEtL27lz559//nn06NFr166FhIR0daEQCET3RyQSzZkzh8fjHTlyZM2aNQBgYWFx4sSJY8eOJSQk4F00ADh58iQeHz8eDwCcnJzwkLt37+Ih3t7egwcPHjVqlFqt3rFjB4vF2r59+/79+9lsNoVCOX369LFjx44dO6bVaouLi83NzadNmzZv3jx8bJNIJE6ePNnf33/u3LmTJk3Cz90eMWIEnoVGo2loaMAwTKVS1dfX451RvHvq7OysVqtzc3PxMmzbtg3DsG+++QYAfH1921gJaM6vy+DxeGfPnsXHQsPCwrKzs3fs2NHVhUIgEN2f2bNnX79+fd++fbiTo6tXr2ZlZZWVleFb4HVuiVo6KmoZwufzAYBGo+3du3fDhg1TpkzBJ/nS0tKCg4P//PPPsLAwPKalpaVAIMAXFctkMmdnZ4FAYGZmFhoaikc4cODAp59+GhkZWVpaamlpKRaLV61adf36dV1fEADwFTR+fn4kEsnExAQPxLfh406UdAV7KWjYs8uwsLBQKpX4oDafz1+7di1aa47oHGoS0squRXZ1KRBdRlJSEgDgA90YhuGKkpWVdfLkyeTkZA6Hg3ez8Gmaurq6tLQ0PCHuugieuj3SzfkpFIqmNgEAX1esVCpPnjyJ9yDz8/N1iz9177qKioqWxcOjffvtt6dPnzY0NPzrr7/u37/f1AkJvltaZ+T1Nk+jnl/nIVOo66Rqc6PGqeOVK1fi3XwXF5fLly/rBtwRiI4m/8i56tgUi0kjurogiK5h1KhRp0+f3r17d11dXU5OTlFREZVKffz48dGjRzdu3Ojj4/P1118HBgb+8ccftbW12dnZCoXCxsZGX18/JSWlqeuijRs36nzyjRo16uzZsxEREYsWLWKxWOfOnQMALy8v3KnQzJkzT58+vW/fvkGDBpmYmERFRQkEAiKRGB4e7ufn5+npKRaLr1+/DgDDhg1zdHQEAI1GAwBEIlEmk8XFxWVlZbVzLbTjODLixZy8VTT964ddXQoEAktY8W1E3xldXQpEl1FXV7dmzRoHBwc6nW5kZDR8+PBbt27hno3x1Z4Yhp07d87X15dKperr67///vsFBQWVlZXz5883NjamUChubm6HDx9uZvbgwYM+Pj4cDofNZnt7e//000/h4eG40KhUqpiYmOnTp1tbW9NoNDMzs/Hjx0dGRn711Ve9e/c2MDBgMplubm6bNm0Si8W4tcrKysmTJ7NYLC6Xu2rVKnz1jYuLCwDMnj0bwzDdCOfdu3cxDPvhhx8AwNHRsY2VgFwadR5nbhdfuCe4tBX5zkZ0MfFLN9fE8sckXe7qgiAQXQaa8+s8iASgkNCsHqLrwVu+XV0KBKIrQeLXeagxUKq0XV0KBAIwRYO2QdXVpUAguhIkfp0HEYBOJXV1KRAIABKBRKN1dSEQiK4EiV/nIZaqquuQyzRE1yPLL5XmFWNaNA6B6Lkg8es8qiUqPRYZ6R+iy6Ea6ptPGCorEHR1QRCILgOJXyeCgZstp7xa8fKYCESHIS0QaFVqmoWJOLW9N04hEP8dkPh1HlEpIgtDWmZRW0/fQSA6AlF0goGfu5G/b9X9uM7JEY2vIt5C0AkvnURqXp2LNdvP2eB6bGXAMMuuLg7i7UKjbIid/bksv7QT8mI5WPf5fRuFw6q6++i233sEYodvv6Fw9fv8toVlh372iLeI7rPJ/c+bxaXCt9cni7JBO3OUpasNZ/PvaXQa6a09xpNIgEUT7Yz0qV1dkJ5F3MKNNh+9y+SZdkJedJ4JmcMCANBqJdlFnZAjkUaJXxbk/9eP5CZu4RD/XSQSybRp03g8no2NzcOHD11dXTkcTnx8/IwZM/CTsv8TdJOe34kbRWqNtq8rt6sL8lzMuTRXGw4AfLvQ43ZCFbytTQ4iEbb8kb5jiSeL0U1+G28/ssIyrbLBdGT/zs6YSGS72nVOVtx+npU3Y3hTR3VOdm8/Dx482LRpk5ubG5lMTklJ6d+/f319fXp6+o4dO/Bjxt5m0tLS7Ozsfvnll6qqquvXr3/22WdOTk6fffYZfibnf4Vu8oKLy6zZt6IXhfzfmMIc09ukq4vwIooqZPdTROP7m3V1QXoKtUnpFH1OV5eiY2Ha8MRJ6Uj8dCQlJU2ZMmX16tWhoaGpqak7d+7EMGzkyJFNfRe8tdBotC+//JJCofD5fBaLhTvnGzJkiI+PT1cX7RXoDuKXUyppUGH/FeV7+7EyYcRl1CLx6zTq0nIYlt28tlm2lhU3orq6FG8RdnZ2uIdYPp/v5eUFAAQCYdasWQYGBl1dtJejE7nU1FQPDw/cscN7773XpYV6ZbqDYOQLZC5W7K4uRffBwYIlEKH9GJ2HRqbQ83Lu6lJ0LCxHazKL2dWleIuYMGECk8kEgJSUFFz8AOA/NGGGk5KSgnv++y/SHcSvVqqyt0D/V+2GpTGdjE5h60TE/Ewyu5v/gGkmhpURD7q6FG8dMpksLy/vP6ofDQ0NmZmZOuX+z9EdxK+4Uk6jdocHeUugUUk5pVKZQtPVBekpUPQ5NFPDri5Fx0LRZ+t5OWtV6DTtZ0hLS9NqtZ6enrqQxMTEixcv1tbWdmGp2kh6erparf6PKjd0D/GjkgmmXHRKb3vSx5UrlqJj2DoJcXIGicno6lJ0OPUZ+eo6aVeX4u2Cz+fb29uz2f/O2tjZ2e3Zs4f+X9gTkpqaamlpyeU2rrGXyWSHDx+OjIysqKjYv38/7pb9baYrxa9SJM7Ma4ddvfnlMgrplR9ErdY8SMx889zfNg6euKZ8Y281ZUKF9LV6fn9culspEr9h7v853rDOWU62ZD3WayfvhDoX1dQfC4l4QyNG/r4qiew1EqrUmv3HQ98w97YQk5ARHZ/RCRnpSE1NbdZzKi8vd3Nzq6mpuXr1qljcef9KFaLabYcuvFKSxMRENzc33cd79+45OTkFBwdfvXp1woQJW7ZskUpfoa1z7GKEsKbulQrwhnSl+FUIa5/kFL+5HRN9GpP+ypNUKrXmflzam+f+tiGXK9/83AInK9bruR5UKBu02rd1D2OH8YZ1XvM4hUh6/VnWTqhzLaaVK990JKAuLQc0r3XOGYbJ5Mo3zL0tqNQalVrdCRnhaLVaPp/fVD8AIDk5GQDCw8Pv379/9OjRziuMGnulX1FkZGRUVFReXl5RUeM5CRMmTAAAIpG4YMECV1dXMpmsVL7Ct6ZQqDr51dEdhj3zy6TEjj+iqUdRWaNUNqDzGDsJq/fHkVjdf9jT6v1xJHb3f8y28+OPP9bW1sbGxspk/3aI8fWTCxcutLa21tPT68LivZgRI0Y8fPjwwoULNjY2usCkpKRx48YRicSKigoajWZo+FbPZHcH8bMyYbBeveeHeAF25ky0hqjT8P7+CxKt+58n5/7NCrqpUVeX4pUpkFVGClM7wvLq1avj4uL27t2L73nA4fP5gYGBAJCYmNi3b9+OyBcA/ii+25Zo9fX1kZGRbTeblJRkbW0NACEhIR999NHzokVGRtbXd/35/m/dCy45RxydWv1KSTbNczM3fMn88J83iyXyHrp8USLX/Hnz1YaXV0xz8LR7Savz0j+C8urOGIxCvJ3ICgW5P59ud7OZO4NVtV3/ZgSAAlnlgsQf7W8t6yDxa4lQKFQoFFZWVhqNJisry9nZWTeo2L7MTzhof2vZCySwvr4+ODh4ypQpp0+39SuWy+V5eXnXrl3btWsXgUB4gfidPn16ypQpwcHBXSuBb5H4JeeI1x1KXfdLam6ppN2Nn7hRNGfb454mgbjszdn2+MSN9v8Xik4RfbQ9bs/ZbCSBzSg+dTVhWVC7m01YFlR86mq7m30NZIWCxOVBt7ymll2NbHfjGTuCb3lN6VoJ1Mne8aI7nZlvUVHR6NGjAQDDMCqVum7dOgqF0kF5FcgqW5VAnewFBwdLJK/wKubz+S4uLnv27Pniiy9WrFjx4rP7JRKJLpeuksC34niz5BzxnzeLk3M7dmmTVKE5caPo4r3S94ZbBgzlkV59llCj0aSkpACAm5tb07XIBQUFtbW1JiYmlpb/Om2pr6/Pzc0lEAgvPe+utLS0qqrKwMDAzs5OF6hQKDIyMgCgV69epCarIbKzs6VSqYWFhZnZiw7Eksg1l+8LLt4rfb1Fm23n5uPKm48rx/Yz/Wisjbnh62w46YjH76Bv6qUUn7qaseOIrKjMaEjvNzTVElmRIGF5UMaOI24bllnPntzu9ttUhkJB5s7gog7WYJVYkrEjOPfn044rZzl8/CF04oRouaZ2QeKPzTTvnujJt5nnOiN7FkCAPZ7XoB8WEMmk4Lp/oCOXQOISGJRx9jPeVFWDIjg4+PTp0001r6ysLDg4uC2m0tLS6urqfvjhBxbrJUuXy8rK8AtcAk+fPj1r1qwGcme4NGlKF4ufoBpbdyi1mezdfFzJz+2oL1wnge8ONldpXq3jK5PJFi9eDABnz551cnLShR86dOjOnTszZ85ct26dLjAlJeWTTz6hUCgPHrzkYIszZ86cPXt21KhRu3fv1gWWlJTged29e5fD+ffU4+3bt6empq5atWr+/PmtWlNriWfuCP6KKmsme+sOtfPQTa7g33XMOgkkKV95LKF9Hx+ng76pFyA4G5a7+zdZUdlrW2gjsqKytktgZmamXC63srIyNjbWBQqFwpKSEgaD0fQMZQzD8HWGzs7Orb68SKKaxOVBzWSvLiUreuKyN3qe56OTQPuPZ5IYrZ9f2I4PiJPUkH9XktIsMF9a0W1cv7VKTYOkUF6p1cgFgspmvb36+vq4uLY6PTY2Nk5PT39ptGZdPYlEIhAImGadfbZ7F4sfm9HK/nQzLs3bsZ2XOTXTVxaDbGJAqyrthgsayUStiQGVxSA3E792r9Ky6uZ7Abv1++G5iPmZjtt/Tql55gcmikr4i9NRqxWgiQT2P7P3BdGCgoKys7M///xzfA0FTkRExJ49e5ydnc+cOaMLbGhowJsLx44da/XMDozBYNhYkPXY6rp/X45kPXa793GF9+OfyReDF6yAb8cHxBnP8PtuxEfHi+4EZZwtlFfhgfNtRgW5zXztJ3oLIfwVgF/ok5mfOk751HGKvLrhV8atbz5fsXTpUnyvHh7BxcWljT2/trN06dKEhAT8evLkyUuXLuXxeL+c6uxN8V0sfnoMmP+u89xxNiduFN2Kq8QDfZz0546zeXHCV0W34sOUS5s7zmZcP1O5oiHtyasZ4XA4rTaCmnZZdAwaNKiNLaZ169Y17YjgODk5tZr8+PHjLzU4prfxZH+LG48rT9woqqxpnJBr9ypNzhHrjL/T13TuOBtzQ9qRM3mvaqfdHx867JtqFX1v17KZk32S02piEnWBel7OvXavfW2brZKyfl9darbuo9GQ3m4blup7u8KTgvbNqFW0TLrbhmWOK2blHjqd+/MZXAKZtjy3De3c88vc+St+QdZjO6780HHFLGAxNL9dad9cXsx8m1HzbUY1k8Buhk72DCgsAJBD4z5OHo8XFBTUTAI7Ap3sdVwWL+atmPMzN6St/7C5BLY7OtnrIPtvG+P6mY7rZ9pMAtsdnex1kP3/BDJH2/5bV0likzN2BIuiEgCAYsAxHtrOPT+KQeO4EC57bbHftOujIzAwsGk/CYdGo7WlBUAx4DSVwDaU+nXQyR7+yCrVczeet/sDNkUnga+U6j/BZtdAney1ik4CX2mrQxsZMWJEUFDQa8ver7/+OmzYsDd3fPhWiB+OTgIrqtvfn866mc49R/aaopPAdrc8uJfR5x+69HDZa4rx0L5DwvsK78dl7GjnYSIdbZe9DkUngc2GKNsF16+W6GTvbWC+TTd0wNvGUVwejzdr1qx2z/0Nbd64cWPp0qUzZ84MCgp6Ewl8i8QPx9yQ1hHv056pfDo64vGnD+uy8Yq3GVwCZYWCdrfc+3AQ0/YtqnOKAcdiyoh2N9vu46iIbsnZs2fPnj07a9asb7/9tumitrZDaK9VTAcPHqyufrXN6UDVAyobJO3/mmgTBBIYOEBN9stj/qcgcF2w2lzAumY7I0HfHpOUgaZn+cLt/nVOogLbEsT5HZjFiyAA1xlqsjo8H4YxAIBc2OEZvYUQ6WBgB9Wdeq53Uwj6DpikFDQvn6A5f/58szWlc+fO3bRp06tKYLv1/B49elRa+mouGjhcU7aeYVlh11Q3iUS285DlpnQ3B5t0ZrxC1mW7g2mMxw0KOYZ1w2W0L6Db1zmBQKTSGUp5lzkkojPjOqGGyRQqAEGt6omHNhCJRI6Bqbi6vKsK0PafcVVV8yVIWq32NXpx7dbzQyAQCASio3n//fcvXryIX8+ePTsoKOj1hj3fujk/BAKBQCBezIcffrh58+Y3WfCCen4IBAKB+M/w22+/DRkypJkfxNcAiR8CgUAgehxvkVcHBAKBQCA6ByR+CAQCgehxIPFDIBAIRI8DiR8CgUAgehxI/BAIBALR40Dih0AgEIgeBxI/BAKBQPQ4kPghEAgEosfRo8UvNKb8nbXRM7993NUFQSAQiHZmZtxewl8Bc+J/6OqCvKW80dmec3fEl4mecaSy92MvHyf9w3/nR6eIHHisrGKJWKIik4lmXNqYPqbvj7DceTIzq0SCpwocZbl4kh2e8OPvkwz1qGQiIadUWiNpYNJIzlbsueNsXK05pyOKb8RWVNU2cDmUkX4m88bb7jmT1cxIg1p7PLzoPl8oqmugkIgWRvRJ/uZTBpnnlkqXf5/UauHH9zdzsWa/yeMjnkdoTPnBi7lGetSzm/t1dVmaI0lKKggKEt+/r6mvp5iams6c6bh3ryg0NHXqVJfg4IoTJyRJSRqJBACGyuUasTh3/fq6R48aSkuBQGC4uFh9+qnZnDkAkLt2rTAkRH/EiJZ3q0JCinbskKalkZhM7jvvOO7dWx8fX7JvX1PLRDq97MgRweFB9X8AAAAgAElEQVTD8txcTK2m8njG775rv307kU6PMjDQiMUtS053dByQk9PJ1dX9yPvyy+Jdu/QGD/aLinpxTDWm2ZN95c/iyFxpOZ1EcWHzjvqt8tazW5t67HLZw4U2Yy6WPXhSV6TCNAO4zg+H7cZTYYBNfLD1emUiAPzqu2Kx7Ts1Kolp+Lwvnaen1BUlivMqlWIOmdHbwCHIbWY/A+cdWSHHiiJK5CIzmsFMqyFb3D+cH/9jfG1unqwcA/jCOeA7j7kFssoNaScjhanChjo9MtNH3+5r1xkjjXvtzw39LPX3Vgv/Z+9P27feuh/tcLC1l70ehUzAr9lMMgBEp4iG9DL6hy9ytGTpMSn5ZdKcUulv1wpYdBI/T+xkyZYq1HVStc5CRY0yp1Q6wYp9K67S3ZbjYceJy6yJy6xNza8b5m18M66SQSMN8jJMzhVfiCwViORpBfXNjPxxvehCZKkBm/LBCKsqsfLm48qDF3MdeCwDFsXNplHhMookAGBqQDPUowAAz5j+Jk+t1mBkEuFNLLyY+3xRyL3SPIFUqwUTLvWjd2xG9zHBWxXj+5vd54sKymUaLeZmw/5xtc/uM9m34pr7ag8YylsxzR61KppS9+BB8qhRWoVCf9gwtp+fSihU5OcDgPDSJXafPg3l5RqJhN2nj/jePTy+sqys4sQJtq+v0ZQp9XFxkvj4jI8+IjGZxtOnCy9f1hs+vOVdSXJyyd69AGA0ebI8J6fq/Pn6x49NZ81qZrk6PDxr+XIgkSxXriTS6aU//1zy/fcUIyObDRv0+vVT19cDgCwtTVNfTzYwYLi6AgDd2rqDqgVTq4v37Kn48095bi6RTme6uLgePcry9o7v04ft40N3cBBevCh98gRTqTgDBvR++FD0998FQUHyvDytTEYxMtLz97fbsoXl5aWuqYkxNTWfP1+ekyNNT1fX1FBNTbljxthv20Y1Ny/csaP82DFlSQnVzMx05kz7LVvS58+vj49X5OUBhll/8YXDd98pCgryN2yojYxUCYUkPT22j4/t118bjByJNTREcblamczt+HGzefMAIGv58rIjRwhk8uCaGhKbDVpttLGxuqbG5fBhi2Xt4A4XA2zG4z1Xyh5xKaw51sOZJBq/rqBKWQcAl8seTrfwf1CTwSbRXdi8J/XFTRNuzwy5UZlIAqIGGr3z/F3+GAAYJFp4RcJAQ5dBhm43KhNvVCZFizLe4/n/UXyXQ6a/a9E/Upi6N+evPGl5THWmn76DWC0VNjR6cQp8vDe2Nru3vsNy+3F3hal3hCmParKqJ560oHMHcJ0BoE4lT5eUAIAXx4ZFpgGACU3vTR5fhakphG7u9qAdHm/TPDdDDkX3MadUWl6tHOJtvGyqPR6CYTB900OJXFNWrTgf1B8AVvyQ1FT8ovgiBo00xd9i3ngbIz0qAGQWS1btT1Y0aG8lVALAygCHcf1M4zJrvwp+Ep1SHfy5n705s6mR0io5AFgY0Yf7GlfWKm8+riQSgc0gWZrQf1ztAwAaLTb+8xgAmDbUYsYISzxVaEyj86obsZUnbhbVSVW9XQwm+Ztv/DUNAC5uGaDHIkenVgcdS6eQCCwGuVai+nC0VVph/ZP8uhXTHKYMMn/z2muVi/cEh//OJxJhsJeRiQGtTKQor1HA01ZFWmE9g0ayMmEUVsjw+HbmTD9nffxapcZS8+sAwMSAiloVzchZvVqrUFh+8onTgQO6QEytFoaGWq9ZY7Nhg+2mTaKrV3USRePx+iQksP38AEArkz20tVUJhdXXr9Pt7RX5+c4//mj92WfN7pb//jsAmM2b53b8uLq29oGFhSI/n2Zh0ScxsalleXY2AJDYbNMPPiDS6ZXnzimlUoqJCQB437qFx0kcMqQuOlpv8OBeV6/iIeKoqKShQwFgsFBINjISXrnyJCCAQKMNKisr3LZNFBqqLCmh29qazZtnvXYtgfLvf+WLwLC0GTOEV66QuVyzOXNITKaEz1dVVSkLCyUJCXbffiv45RcSm81wcZE9eaIrPIFCMZo4EVOpqq9fF16+XBcb619SIvr7bwBQCgSS5GQ9f38CmVwdFlZ+/Hh9QoLBsGGlP/1E4nCM3323NjKyeO9eeV5eXUwM289PIxarhI3+Y9MCA+tjY9m9e/OWL6+9e7f2zp26R4+GVFcTqFQ9f//aiIiau3dx8auNjMS/O3FUlOH48ZLkZHVNDQAYjBjxOr+MFoRXJFwpe8Qh05NG/mDDMNGFJ4rz8mWV03kDBxnOB4B1qcebit9dYUpQxtkvnd8LLrgpUjVK1yXBg5HGvaaY95tvM4pHNwSAx7XZ/e+tl6gVJ0siAeBgryXzbUbdqEwc/2DLpbJHKSMPeOnZ9I1cpxO/bKkAAAZwXWbwBmMAkcJUQyqbRCAGWg4JtBwCALerkt+JCQKAY73/19eg0b/PsaIIAFBjmrWpx44VRZAJpAU2o6eY9xsatQEAhBNOGFE5V8oeBcR+RyOSDSisCqV4g8t7D6ozo0UZ+3st+th+fLvU5FtLO4jfol3xKjVmbkifMMAsYBgvii8y0qN62HIA4MbjyuJKWXaJVCLXmHJpEwe0rhZRKaIB7lxna5YuRKHUAAAGGGgJAIBb87Dj4HfT8uvszZlNLUwfZpmYLU4vrF+6NxEA2AzSprlutmbPxHke1fUNv10rcLVhP86oiUmtplEa50EJrb2Bz0SUuFizh3ob6TE7qllUJ1X/HlYIAJvmug3pZaQL17UqPO04AHDk73yd+H0w0vKDkY2KfvVBeWp+HYtOmjjQ/PqjCtSq0KEWieofPwYAeVbWAwsLjVyu16+fw549apFILRIZT5/eMgnF1JRiaopfYxqNVqEAAKarq/DSJSqPZzhxou5Xorurrq4GAH1/fwAgGxgwPTwkCQni6GjeypVNLZsGBpYeOiTPzEwcMgQAgEh02L3bYsmStj5M018nhiUNHy5NSWH37m3z5ZeV58/nf/WVPDPT9dixtliqDg8XXrlC4nD6JiXRbGx04SU//EDicAzfecdo8mQAyF23Tid+VmvXWq1d2xht//7czz5rEAi0MlnVpUsGI0fafPUVp3dvIpMJAMV79uStXy/l82VpaQDgdPCg+fz5NTdu8MePF1661DclheXlFd+3r0788DaB3oABJjNmAIbVRkZSDA0JJBIAcEeOrI2IqL17FwBUFRXyzExWr17SlJTayEjD8eNxLaSamzNcXeWZmfmbNtU9eKCurWW6udls2KAoKCg5cEAlFJoEBJD02tQlCq9IAABjql7g4738ugIzmsECm9EbXN6/JHjIo3P9DVvxpFOurPkwbt8wY4+t7rOCC27igVKN4mZl0v5ei3rp2epiStWN/nI1GAYAuLVBho1uCmKqM7z0bJpa3ugyY/2TP34puP5LwXUA8OLYXPX/mkRo04qNkNIHzmyLXnq2/4jSdudcNqY2Pn6r77edWRf7Gji9x/M3onLaYvw/zRu9wRk0Uh8XAyN9alGFLKNIcvjvfJlScz9FOMjLEK/ZiPjKxGwxABAI0M+Ny9VrpSlaU696UlC38aN/f0yFFbKdp7IAYKCH4aO0GjwjAGBQSXiEerm6mRF7C6a/l2FStvijsdbFVfKL9wS7TmcdWuOLv/FfDIbBvhW9bMwY+0Nyrz0oL66U4+GE1n4dAzy42xZ5vNTmm5CcK25Qa4kEuPm4cu/ZbAqZ2N+du2SyXdNWxfPQaLFzd0oAYNpQHotOQq2KpjSUN0qyOCrKfP78+oSEmtu3+WPHGk2ZwvTwYL7QQ4q6tvZJQIBGItEfNsxy1aq43r2Np03TPYnuLmfAgPpHjwCAxG7sFpM4HADAOyVNoRgbm8+bV7h9u+2GDSQWK//rr/M3bmT7+nLfeadND0P898WHabXSlBQgEt1PnaIYG+v5+/PHji3/4w+H776jmJm91FJ1eDhenrTAQAmfTzUzM1+wwHbDBuGlS0aTJhFotFZTKYuKBMHB6pqaqgsXgEi0/vxzwLCamzed9u/Xx+UcAAA0UikAEDkcbX09PG0T6A0ahN+ti4lheXk1NWu7cWPu+vWCX34R/PILALC8vHpdvQokEjzt0imLihR5eXgjxmr16tx163DZw/8ajBghz86O79tXI5EYBwSwPD1Lf/75ydNmjeGkSQ3l5TVnzrShfqFcWQMA+bJKK4bRApvRx4oivsk4owXsUtmDaRYDCND8R6zBtDMf7yMRiGf6rm0qS2EV8Q1a9TSLAbqQtPpifBHKZPO+VyviAIBNpuv+AkCNStLM+FhT3wuCaBKBuNRu3IXS6GsV8QsSDkYM3tKyGC2xYhgljfyBRqS4RazMlAgyJCV4OLG11Y6TzPqGDtz4Upvdgzd60fyyxpf4tPI3H0uPSa2+/qi8srZhVYAjHrh7uZdKg+UJpJt/T7/2oBzDsM9mNHe5G50qIpMIA9wN8Y9RKaJdp7OUKu2SyXY2Zgxc/ORKDQDIGzR4HA6jebF3nMyMy6xdNtV+kr85ADxMqy6tUtyOqwocZfnSp9BnkW3MGADAZVMAQPE0FxyN5hmXT31duC81+IZU1zcAgBaD/DLpuP5mUSmim48rK6qVNZIGXavieUTEV5VXK+lU4nvDeKhV0QydIPFWrHDYtUtdWxvN5aqqqoRXrliuWPGChLK0tNRp0+TZ2SYzZridOKEoLJSlpTkfPNjyrvPPP8eYmgIAvrAFADT19QBA5jb/2RTv3Zu/YYPRlCk2GzYAgDQtrSw4uOSHH9oqfgAAgKnVAEAAwABAq33s7t7kHibPyWmL+OFtAkV+Ps3KynzBgvJjxwq++UZTVyeOifE4e/Z5qRRFRUXbt+PXDBcXg+HDRWFh2oYG42nTdHHwhT8ECoW3eHHJDz/A069A90W0bBNwx47lXLhAIJF4S5dWXrhQfe1axoIFPhERQCBw+vcnslhaqbT27t36+HgAMBg1Sn/o0OrwcE1dnfj+fQAwGDFCcPiwRiJhuru7BAcDAJBIhd9+CwDG773nGRICAEkjRugGn18Am8wAACIQwgZ+wybTKUTy/tzQk0WRObLyg71a6aDnSMvuiZ64sHiLEn8CgDq1DAB+zLum0qoHGbqZ0QzwaJfLHs6NPyDTKHd7znVnW+PiJ1ErdH8BgEt5ZspcqlGMjNokUtXHj9jbW99xrIkv78bCu8LUxzU5/bnOL32Q/lxnGpECAGY0g0yJQPa004mjxp553Y019X2pwW7D64tfnVSt1mJNZ/sAQKLQcJhkb0c9qULDoJKIRKCQCK7WbJ4xXVTXkF8ma2knii/q58qlU4kYBn/cKDp1q5hFJ21b5NHfnVsrURGJoNVCemG9tSkjraBxBNzDrvnARa5ACgBKlQYAMAxUKgye08NoCZn0TI+E+DSZRK5mM0g5pc+0wqiUf5tL9/kipUrrYcfhGb3RLFczdFL0v+mO/d25Po76m4+lJ+XWEoCga1W0ihaDMxElADBlkAWHSb76oBy1KppCs7WlWlg0lJURiEQAIJBIQCAAhqmrq43fe+95qYSXL2fMm6eRSu23b8eFquriRbKRkf7w4a3excfi6h4+tFi2TFNXhw/36fo6OiRJSQCglTc2CLQyGcDLf7JNZYNsYCBJTARc+QAIZLLnxYv4YCMAiO/ckaSmUs3N6Y4v+s38a5NI7BUWRmKziRRKyf79FadOEalUo4kTn5dKf8iQ4RimEgqLd+0q3rs39d13jcaP1x80qFFutdqCzZsLt28n6+t7hoRo5HJc/PA2ga5l0KxNoJVKk0aOVItEfeLj2b17c8eOfcDj1d69W//4Mad/fwKFoj9oUM2tW7WRkfXx8TRbW7q9vcHIkaLQ0NKff9ZN+NXcvg0AsvT0GBOTpsY5fn66i7aI30Cuy/GiOwCAd+NIQAQApVZlROEMN/ZsGR8DDACypIIsqUAXyBcXkojEvXbj8AibM85sy7ygT2Fd7ff1BLPelUoxmUBUY9qHNVmubMsH1Zl4qmZjqiVyET59KNc0AIBM06hebXy/UYmN/9R4N5Hy9GONSmJAYSWK85pGZpD+bdReFDyQaxr8DV0dWR21uKFreX3xK6qUrT2U4mWnZ2ZIL6mSpxfWAwCDQurrZkAiEqJSRL9dLfCw5XCY5IJyWWaxBAD6u3O/P5+j1mDl1UoAeJRWU1HdkJgtXjfTCQD+ii47dasYAAw4lJB7pSH3SgHA18kgIav20JW8xxk1STliABjsZXj5vqCpkeo6lbEeraZedf5OqUyhEQgVlbVKMonQz/113qdkMoFFJ0kVmr1ns+3Mmdcelj8v5sGLubUS1WcznNpX/NyfDmwSiQR4Or5FAALeqnhBwn+ShSVVciqZiM+9oVZFS2y//jp75crSX34BAqEuNhYwjGJoSNLTY/v6AkD1tWuVFy4oSxrHhbKWLFHX1YlCQwHDKCYmNXfu1Ny5AwCyJ0+Mp04lkMnSlJQn773X7C7DwUGaklL+xx/q2lpZdrZWoaDZ2tLMzTPmz29qWVFYCAA1ERGZixaRWKzKc+cAoGm3qVWYLi5kAwN1bW3mwoUsLy9BcDAAEIhEppeXNDW1YMsW46lTtQqFNC1NdPUqaDQuv/5q8TLx0xs4sPz4ccBbAwD4GKNGJuOOG0dksVpNohaJyEZGAEAxNjacOLF4715MpaqOiLDfvh0ANGJx2uzZ1deuMT09va5cYTg5qSorCWQyplbXPXzIcHWte/CgMWt//6ZmlSUlapEInrYJGhsE8O8PzmDkyJpbt0RhYerqanzZCz4WWrxvHzyd8MOnLZnu7k5Pu+a5a9dK+XxFQQH+UZ6fj18IL17UyOX6/v6ttg9mWw3flX0pX1YZELvTV98+uPAGAJCIxKmm/cgEEgDsy/krpa4wrjYHAHKl5d9lXVpmN/awz8d4cuOwuSJV/f8cJ/6YFxbAGwgAP+eFb828AACmVP19OX/ty/kLAEaa9LpVmbya/9v1ioQ7VSkAEGAx4GDeVZVWUyCrBIBr5fECeY0BhVmrki1M+CnQanBYeQIAODDNmk4ith0uhWVAYdWqpAsTf/Li2OjmJluykn+kQin+1XcFEr/mWBjRx/c3S84VZxZLyCSCizV7uK/Rr6GFQ7yNAMDGlGFtykjJr5Mq1BwG2dOOM7af2YQBZpO+iFE9bfIXlMvyy6QEImGgpyE8nY4CgNIqRWlV4wjAggk2XvZ6N2Ir/uELuWzqe8N5CybYBmx82NRIQbmsj6vBe8N50Smiy/cFdCrJ20Fv9jvWzaav2giJSPhytsvhv/KzSiQYhgWOssIludOwNmUM8zH6J1n0y5W8oT7GkUlVAKDHpPh7cklEAgBciCzNL5NlFUsAQCBS7D6TTaUQV7/nePp2MQBMHGjG5VAkck1SjnhtIGpVPANvxQoClVryww/F+/ZRTEzM5s6tvXdP1+2TpKRU/PGHLnLFyZO6a1VVVW1EhO4jvjpGI5UChjW7yx03zuPcuaKdO0VhYUQGw+T99x337q04c6alZat162pu3aq6cAEwjOnuzlu69KULXohMptuff+auWVMfH49hmM0XXxRu3QoEgu+9e0U7dwr//rto504ig8FwciIxGLoO1osxnT27aNcuRX5+akAA29e3LDgYALQSicnTqbKSffskKSn1cXEAoMjNzZg/X3j5MsvTk+HkpFUqq2/cAAASh6OprzcJCACAjEWLqq9dAwASm521fDluxGTGjMozZ7JXr66+fh1vJRgHBJQePKhVqXBZqr52rUEgwKU9Y+FC08BAUVgYANAdHFi9euFGcKnDlxTh12wfHzKXi0smHsJbvrwsOFiWni44dIjTt69KKMQXIpWfOEExNtYqFKLQUNxa1sqVqoqK57UP2GT6vSHb16Uev1WVFCVKd2Hz5lmP+jT16HSLRsG+VZV0o7Jxt4+wof6P4rssEk0nfjgJtXl99B1tGSYAINE09vKb9g63uc8abOh+vOhOiCDGlKa/xnHqNvfZ3LDZSm3jBERqfVFqfdFwI09zusF9UdqurMvGNE6g5eCt7rPoxJdPPbSEQiT/2fvTNam/x9fmYBj2hUsALsk9EAKGYS+P1TYu/SM4fr3o4pYBum1/L2XzsXRFg3bXslaGEXoy+O66u4lVYonKlEvr5869cr9s6yL3gR6GAPBV8JO4zNqm8elU4ldzXDf/nk4mEU5s6GtiQL0VV7n3XHbIlgEcBvlsRMnRsMJmWSyYYKPRwo3YiiqxksumjvAzbtaqwOnjamBnzoxOEYnqGuhUkr05c/Y71r1dDHQRdKs9l06xa7ras+km9xM3iv68WexsxZo7zubwX/nV9SpnS1YvR/1Tt4p1qz0/m+E0cWDjBNWMzbHNQjoIaXJynK+vX3R0y2HJ51F64ED+pk2Dq6qetxLkv4iyuDh33bqaW7e0DQ1MFxdO//7lR48OqqzEhyX548fX3LjRND6BTKbb2SkFAgCgmptzR49WV1fLCwr6xMW1Gh8AfO/fr42IKD9+XFlSgp8tYL9tWxSXiymfmYLSHz6cam4uvn9fJRRSjI31hw6137qV4dw4s4Wp1dFcLi7qA/Lz6XZ2AJA6bZror78AQLfDT5aRURAUVBcT01BRQeZy2d7eVAuL2shIjUTCHTOGxGKVHzumN3iwPCenUfwWL25LLR3Iu7op/VTVhBP4FNpL0WBa8+vzP3OcusHl/bbER3Qy7Sl+95KEDWrtO31N257k3J1SdzuOt8Mb7cfs9qBWRQdRHxtbffOm7caNbR3JBag6f16rUJjNnduhBetangQEaKRS75vPHRBrjkYTY25u9dln+Kxnd+V8abRC2zDXemQb41cqxYfywxfajm66TRDx9tCe4ofoIFCrAtGZFO/erTdwoP6wYW2Mr6qsLD10yGLhwqbbBBGItxwkfggEAoHocfRorw4IBAKB6Jkg8UMgEAhEj6NDjpK6+qD8aky5QKTQaDAjfeogL8OFE22p5H+FNlcg/e1qwZOCeo0Gc7FmL5pk62Gr18zJwIIJtk3Xd5RXK38PK0jOEYulaiad5MhjzX7H2tdJXxehpX+lD0db4Zu+dXg76O1b2asjHhmBQPQQnufqSBchua7gyycnoqvT1VptX67jDveP/A1dm7ku2uY+W7dq9HjRnQWJPzbLxUfPLmnkv674nG59nCt7ZnfQBuf3d2SHNA0ZbuQZOWRb+z5sN6b9xS82veZASC6RAFOHWFDJxL+jyy7eE+gxKbPGWOERyquVa35KkSk1LtZsBpWUnCv+/JfUwV5G95KFTZ0MlFcrvpn373GL205kZBZLnCxZkwdZJOeIk3LEGUX1l7YNpDzrAaCpfyUalQQARnpU/JwRAHDgtb5p978CalW0I1qFomDTpqpLlxoEAgKVynBwsFi2jPd0XxpOS898RAYj74svRKGharGY4ehotWZNs4Xyz/PSh99VlpQ8bOGZyHzRovKjR5uG2G7ebBcU1P7P3Gaauxby9rbduNFg9GhdBJVQ+Ew9fPqpxdKl0uTkvC+/FEdHY2o1p08f++3bm66aeabCKRS6gwNv6VJek4PlKk6ezPjoo6bF4AwYoBIKFbm5TQN7hYcbju8yhwMvcHWEky+rGHZ/Q51a3s/AiUWmRwpTR0d/E8AbcK40uqnrogJZ5YV+6/EkPLrhaGNvnYUHNRkyTYM1w7hl7kON3KmERsnEXRdZ0g3d2I2vVh99u4556O5J+4tfqVAOAAwaabiPMZVCjEwSKhqUtZKGd9ZGA8D1PYMu3iuVKTXmhrQfP/EhEuGTg/z0wvp7fCE862TgPl+UXy7TbVTHzbrbcob7GANAcq6YzSBNXB8DALuWeep2njX1r3Q7vgoABnoafvp+K5tYD17MjU2vqZWo1BrMxIA6ws/ERJ/246VcIz3qzFFWp24Xq9TaaUN4Y/ub7jubnVEksTJlrPnAybXrPNWhVkX7UrB5c/HevRQzM+t165QlJeXHj2d//DHb21u3568qJCRtxgxo4pmvLjaWoq8vSU6m2doaDRkiCg3NWrIEUyp1Hhte4KWvadZEJrPZ4SYAwPLy0h3FyXBw6NiHfxmNroV8fS0WLxZHRdXevVv38OHg6mpcxTG1mj92rCQxkWZjYzR0qCg0NGvZMmVJScmBA5q6Ok6/fiQ2u/bu3eQxY/yiozn9Gh0a6yrc8n//UwmF5ceOZa9cyfLyaraslMrjMZ8eUsry8MC3uusPGaLbWEkxbkUVOo1WXR19mXZiTMzmwYZuUUN3/pD7d51a7sA0ezBsF4lA9P/ni4fVWedLo+FZ10UhggepdUW494axpr66QzXT60u87nwCAJ84TiL8FQAAtwYFjTHxwe+e7/e5Oa3xiImTxfcAYLJ532Y763FWJB8Jr0yoVNaqtBorhtFMy6FWDKOV/GBLuuEXztO3Z11QatWr7CfOtxm1KPGn2JpsFzbvV78V/Qxeflhot6H9xW+Er8nf0eUlVfLPfkoBAAIBlky2s7dgXr5fhkd4UlAHAG42HPzgLg87TlphHQF7iZOBWWOsf71aEBpTjrvLsTNnfjPfbeF3Cc1yb+pfSY9FAYCI+MqbsRUcJtndlrNgoq3OI0GuQOppr2dlwqiTqcIeVJyNKMHPpqmubzgfWeLAYyVk1Z66XfxXtMDegmWkT80TSHeeyjz+ZZ92r7E2gloV7QvuPYdhb2/ywQfK4uLy48cJZHLWxx9L+XybjRvtt20r3LIFnvXMpywowLdk+969S7e3z//qq6LvvivYsoX38cf4MXSteulrqKi4RyAAwDCVCs+aZm3tc/u2riSZixcDgM3GjaYzZzYrpFYm40+cKM/MVFVXE4hEppub1Zo1omvXqs6dM5kxg2xgUHH6NMXY2GnfPhKbnb16tbKkxGD4cLdjx3SemN6kcjgDBpjMmEGgUGrv3iVxufcZDADoFR6ulcnwY0V9IiIYTk75X39dtH178Z49WoWCbmfnFxNDIH9zyGsAAA8/SURBVJOThg4VR0UV7djheflyywpXVVWVHzsGJFL2qlXSlBTcky0ejWpuznRzoxgZ6Q8dyh0zBhc/lrc3gUxmODmZzZ5NNmw8rrbs119L9u9XlpRopFKKkRF3zBin779/7OOjqqiw3bRJHBNTFx3N9PBwO3685tat4j17tEqlaWCg808/4Ue4vR6tujrC4N818zGiTADoz3XGDwUdZOj2sDoLv/1S10UAsCMrRAvYUCP3oUatnPPuEfE/pVZlzzRbZDsGd050qvif40V3DCmcgYYu29xne3AaxxWS6/IHG7q5sHmihvpfC27tzL74Hm8gAJQpavbkXPbWs7tVlbwt68LP+WG99Gx5dMPkuoI5cfszx/z82jXzn6P9xU+fRRnbz/T07eIPR1vTqcRj4YW/hxWun+W8bqYzABAJBIlcA0/9CQAAk0aCp445mjoZwDA4E1GCOyj/7AOnvq7cf5KFRCJh0kDzf5KFj9Jr9l/IXRvoTCCArTkTWvOvNKSXkY0pw57HIhIgNr0mOrX6SUH9r5/7GbApALDnYy8qmYhhoFRpNRosNKa8qFyG5/v9Sm9zQ9qsrY+rahucrdi7l3sl5Yg//yW1tEohkWvYjNf/z3kTUKuifbFavbomIqLu4cM4b28AIHO5HufP53/9NX5XI5FIU1KghWc+AKCam9Pt7eHpcdWqykp5Tg7DxQWe46WP1atX6dNzJnEUeXlRHA6BTGZ6etp+9RUemLV8ecb8+TQez2jKFLvNm/FXvLahQZ6TYxwQQLO0lKSkVJ07lzF/Pnf0aACouniR7evLdHaWJCWlz5lDoFAMRoxoKC+vDgvL37TJ5ciRN6kc3LVQ2ZEjZUeOAADT3d3z0iWd1whxTAwAUMzMGE5OAKA/eDAA4KeIcQYMIJDJeOWIo6LE0dEtKzyhXz+8Sj3Ons3fvLlZ1pKEBLyeAcB80SL8QnDoEH5RsHmz7717+IFnsvR0uo2N8bvvEsjk8hMnKk+fVovFeLTCHTu4o0eTjYwkCQlJw4cDkag3cGB1eLjg8GFO//7mCxa8ds206upoqe3YY37/w1031Kqk0MQ/EYfM0KV9xnURBjuyQ04U3wWAYN8VuGjlycrPlNwHgK9dPqAQyMf8/gcAnhwbPNVYE18ewzC9vvhRTfaa1GPTeQPd2Vbe+rZEIIZVxF8uexQtykgZdcCUpg8AEYO30IlUDDCZRqnGNIfyr6fXlwAABtg/Q3bYMU1tbi4ulov6GDjeGvTtXWHKqOhvsqSCWpXUgPLfG8V5Pdpf/M5HlvweVjjQwxAfjiuskIU9rLgdV7VjSWNDBhcP+dOTPGVKDQCG698zTgYIWEWNsqJGCQA19aotf2TUy9SHPvN1tmL1cTWY+e1jfq54yWQ7nT/xlv6VckulJzY2vlJzBdLl+5JqJar4zNrRfUwwDE7fKrkVVykUK7VPG20ypQYADNgUc0MaAOizKFW1DW42HGji7kCu7DLxQ62K9oXVq5fx1Kk1d+7Ybd4sy8ws+f77jLlzHXft0i5fzvbxUdc2HiDXzDNfqyFPZswgcThEGs3n5s2WXvrc/vgD9yuLnxxNNTfXGziQpK8vvn+/Ljo6ZepUo4kT9fz9Gc7OKpGo+vr10oMHpXy+z927AEDW0xuYn0+gUDC1WqtQyNLSpCkpDRUVAMBwdOwTG6soLHzk6KhVKOy//tpm48aCzZsLt2ypj419w8pp6Vooa/FitxMngEBg+/hUhYS0Wg8tA1VVVXg7QK9fP9tNm4ynTq2JiLALCsIrPH3ePKe9e7FVq3Cvfiwvr17XrukNHAhEYkFQUOmBA+VHj9ps2GA+dy7N1lbK56dOm9ZQVpa7di1++oz9jh34MKxWJmN5e6fNmFFz+zbZwAAAzOfOdf399/LjxzMXLFDX1PhGRuoPH85/552a27frY2PfRPxadXUUXZ1+xLdx7BEXD51/onq1XJe2meuiQllVoawKnvo/AoCdWRc1oO1n4ISPgs63GaVLmzByn84DX0Dsd1fKHiWL83PGHMZDkusKfO9+VtkgvlWVNNtqOAbY9syQE8V3S+XVGtDicerUcgAwoerbMU3xi2K5CHeKpHOiJFErkPi9Prml//oBAAClSgsAKrUWHysb3dvE004vu0SaVSzRaoFIhLSCegLubQN7xskAAQjB6/zsLZgAUFwpr5epdWZxmxjA44yakip5bxcDMpHQ0r+SvEGNYc0PrmpQawEg5kn1qdvFZBJhTaAzz4j+d3RZZJIQF0HSUwnFE5JJbT34qqNBrYr2JX327Orr1x337bNYuhQARFevyrOylGVl1uvXw1M/fNDCM1+rIVI+HwAINFqrXvoq/vyzV1gYHpNmaekvEOC/La1M9sjZuUEgoJqZeT09bbn0wIGcTz+tjYxsKC2lWlqqxeLsVauqr1/HT3PGwftYbD8/IJEoTx33cPr3h6dOgtRPC/Z66FwL9Y6N5fTrZzhhQoyZmTg62mH3bryziwtMy3poNbAuOhoAiHR62qxZNTduOOzejVd4dXi4LD1dUVSkmxBl+/qCb+PUl/3WraUHDgAAjcdjuLriD2j64Ycl33+Pj7gCgPDy5aLvvpNlZGANDXgIplTiR43jE40UI6N2r5xWXR0RgHCy+J4pTX+sqe8gI9d4cW5cbY4G05IIRNxREe5z8RnXRQTgj9zf1DNDiVx0oigSAL52nQEAWtCeLr4PAGNMvSkEsgpT62b7cCQqBQZYM3+2Co0KAP4qi92WdYFKIP/mt9KRZX4oP/xsaRR+ngmF2Phf1ujhiNBRvqPfftr/yf2cDSKThInZ4r3nshlUUmSiEADsLZi7TmcBwEg/4+nDeNdjKwQixac/8elUYnphPYlIGORleJ8vauZkAFc+ALA0ZpgYUKtqG/aeyx7haxKbXg0A5lzqiRtFALBrmSeVQmzpX4lOJc3ZHufy//buNabJM4oD+OnbC73QlrVSKdoiggTXYmsm07XOqTgNW5ibkbilc8mWbbhlGjVmGtsZjbtkTqbZEpNlEshgxBm3EUNWUYPzFsJtjoAQCrZQK9AC5VJqS2nfdx8eLNjxSc2yref3DUKfwCHhnPI8z/tXJXLZrIaOEQCQiDjPaWQAQNMMALAAwmHa6b7f/PBzov+dcKp4sh5K1GMY0lECNpu7oiIxJ0ek082ZzAcAoYGBYE8Pf9EiEs3DVSgM/f1kz6/jjTfgbyl99OQkiXGYbzIFHQ6eUkkJBAAzST3knVwM8v30HDniqawUZGZml5ZyZLLOd98NdE6nvrG43NmLTH/4JMREC0UeRAt5a2oCdvtTeXlSg8FVXDzldgftdv7ixaQOFJ9PB4O+hgaIRIDNJp+Ub96sraoiL69TKuHvBe/qcldUiLTaRL1+tLZWajSSgy3ktwMA4fGZg5QTLS0AwJZKASDs9Xa8+SbQtGr/fnl+vq+p6c6+fdGvnLM4DE0/fnHmjDpSC+dt/+OkUZa9UaHfnVFQ2lvb7R94/vpBESehbqSTS7E3K58911cXE10Uk0n0ZdcvISa8TJJWkJILACE6vP2PkwBwyXCYT/HW3bQYZUsXCRW2ib66kU4AEHL46ReLViRl8iiO1UN2IsXktRGGJj99iA53+FwXH6RPoNmefPN7adX80BRtbXBfaxlmGEY9X/Dycympcn50a0op55/4aNnp6p7bPeM0Ddp0ydv5aZp0ceVlV0x0UXRNioJjO7RlF5xt9vGfal1SEfcF/TzTBtX7x29F14zJVyowpFAs1sUmT5t9fCIQlkl4xhy56UUV+decMUe+ZU3qpSbPqSqHfol09TK5tX6OP0D/KjhVPFlJ69d7Kiudx46Fx8cD3d2TTieLxxtvbOwvKVGbzek6XZrF0r5t20PJfGo1Ryr1t7a25OWJV6wY+vVXAEgzm6dDF8maZ87EpPSJtFpygl/x+uvuykpXcbHEYOAmJ4/duBHq6wOK8lqtzcuXCzWayNiY98IFAJCuWTMdMheJAABQVOT+fV9TU8Bm+wcqI8jMTFCpJu/e7XzvvWi0EE+tJieAcqxWeUGBSKfzt7S0bNggzs0dqqoCgIV799775pugw3Fr9Wq2WDx2/TqLw4nuaEYLfre4OOL3B7q7gw4Hi8v1NTcPlJWp9u9P1Ot7P//c39oqNRhYCQnD1dUAwFUoHGbz8Pnzgqwsf1sb2QtcuHMnPEixJyUKud3uH3/8ByoDc0Ud7c3Y3Obr/c09vU+5WJhy7fnPDtwuv+ntoBlmtWzpp0+bjLJsre1cTHTR7GXdk6MlvZcBwJxVGPNmDgAyRClvq/OuDt1uHO3isji5SZk70jdxWOwyZ+0Nb/tIyK/kP/WqcqUlq5Bs+L2WunJ3RkG58/fdrSXrkrVbUled7r0M6GH4bM//kqrr/dYGd/9wkGEYpZxPpoqD37cDwIWvDGyK1X3PH50qslSJs6eK2dFFCbOSY12DATJVjPmnpCKudrEkOlV8WaRJSxH+UONsuTM2NBrisFkLkgXRqaJ34D6ZKvSZUtOLKpK9RzPw3XnHpSbPZIjWL5HKJTxrvVsu5Q2PhaIJRx+e+LPL5d++UfXWJvWde/4dX/8JAJWf5CYnPUo+2SOL+Hw9hw8PVVWF+vookUik1aZZLA6LxVdfT057AsDg2bPOL77wt7dTAoGM3PMTCu0ffzxcXT1zv62oaPay9779tr+kJGi3A8OQq2z8jIzW/HwAWDM15WtsvHv8uK+xMeTxcJKSxMuXqw8c8NbUeGtqgnY7HQrx1erkwkLVvn1siQQApgYHO995Z+TKFYrHm28y+Zqbx+vqBFlZAZtNYTItraiITEzcEIsBQHflStLata6TJ+/s2cPPyFjZ3f04xQnYbI5Dh2aihYzGNIulSaeDB9fspgYHZ9dhwa5dqR98MHHrFrnnB5FI4jPPpB89mrRuJgMhtuAaTZrZ7Dh0yFdfT057usvLB0pL/e3t4dHRhAULZBs3KouK+k6dGr16ddLlohIShNnZC/fsSS4sJAsOlJX1HjkSGhwUaTTzXnmFHFbiKhRTHg8JKhqurm4rKCCVZ3E4t7duHfr5Z1K3xykO+n/A5ocQQiju4LM9EUIIxR1sfgghhOIONj+EEEJxB5sfQgihuIPNDyGEUNzB5ocQQijuYPNDCCEUd7D5IYQQijvY/BBCCMUdbH4IIYTiDjY/hBBCcQebH0IIobiDzQ8hhFDcweaHEEIo7vwFASTSc99QpdoAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.011669,
     "end_time": "2022-01-02T18:59:44.351576",
     "exception": false,
     "start_time": "2022-01-02T18:59:44.339907",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# The sampling strategy\n",
    "\n",
    "The idea is to sample from the past segments of the data that are predictors for the future segment. \n",
    "\n",
    "In the case of Traffic data, \n",
    "![create_data.png](attachment:66f43f1d-1b32-4eb8-8c59-814f64b9d4ae.png)\n",
    "\n",
    "If we want to predict the traffic in the next hour e.g. (Thursday 8:00 AM) \n",
    "\n",
    "We can use :\n",
    "## 1. The recent segment\n",
    "\n",
    "We use the last two hours. Based on the traffic from 6:00 to 8:00 AM, we will predict the traffic of 8:00 AM\n",
    "\n",
    "Intuition : The formation and dispersion of traffic congestions are gradual. So, the just past traffic flows inevitably have influence on the future traffic flows\n",
    "\n",
    "## 2. The daily-periodic segment\n",
    "\n",
    "We use the same hour in the last week and the same hour yesterday and the day before. Based on the traffic at 8:00 AM on tuesday and wednesday, we will predict the traffic of 8:00 AM on thursday.\n",
    "\n",
    "Intuition: Due to the regular daily routine of people, traffic data may show repeated patterns, such as the daily morning peaks. The purpose of the daily-period component is to model the daily periodicity of traffic data.\n",
    "\n",
    "## 3. The weekly-periodic segment\n",
    "\n",
    "We use the same hour in the last week and the same hour in the week before the last week. Based on the traffic at 8:00 AM on the thursday in the past two weeks, we will predict the traffic of 8:00 AM\n",
    "\n",
    "Intuition: Usually, the traffic patterns on Monday have a certain similarity with the traffic patterns on Mondays in history, but may be greatly different from those on weekends. Thus, the weekly-period component is designed to capture the weekly periodic features in traffic data.\n",
    "\n",
    "\n",
    "\n",
    "### The three ways can be used to build three different models and the output prediction of the three models can be fused as the final output. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.011628,
     "end_time": "2022-01-02T18:59:44.375398",
     "exception": false,
     "start_time": "2022-01-02T18:59:44.363770",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "\n",
    "\n",
    "In our example we sample data per hour, so ignore anything about weeks or days\n",
    "\n",
    "Each hour's data is used to predict the next hour's data\n",
    "\n",
    "the 1 hour data(INPUT) is of shape (12, 307, 3) and the next hour(TARGET) is the same  (12, 307, 3)\n",
    "\n",
    "The functions get_sample_indices and search_data simply finds the input and the target since the orininal data have 16992 points collected every 5 minutes \n",
    "\n",
    "\n",
    "for example, \n",
    "the first 12 points(0:12) will be extracted from original data have 16992 points i.e. original[0:12]\n",
    "The target is the next hour so it will  be the next 12 points(12:24) i.g. original[12:24] \n",
    "\n",
    "Cool ... generate more .. move 5 minutes in the future ....\n",
    "\n",
    "the first 12 points(1:13) will be extracted from original data have 16992 points i.e. original[1:13]\n",
    "The target is the next hour so it will  be the next 12 points(13:25) i.g. original[13:25] \n",
    "\n",
    "\n",
    "And sooooo on. \n",
    "\n",
    "How many INPUT data and TARGET examples do we have now ? \n",
    "We havee 16969 examples which is 16992 (original) - 23 \n",
    "\n",
    "The missing 23 points are the first 11 points(first hour) in the original sequence and the last 12 points(last hour) in the sequence i.e. the first hour data can't be predicted from any other previous data and the last hour data can't be used for predicting any future data because it's the last one in the training. So we ingore both. \n",
    "\n",
    "## Data splitting\n",
    "\n",
    "10181 data/target examples will be used as the training set ( 35 days )\n",
    "\n",
    "3394 data/target examples will be used as the validation set (12 days)\n",
    "\n",
    "3394 data/target examples will be used as the testing set (12 days)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:44.410034Z",
     "iopub.status.busy": "2022-01-02T18:59:44.408973Z",
     "iopub.status.idle": "2022-01-02T18:59:44.411944Z",
     "shell.execute_reply": "2022-01-02T18:59:44.411362Z",
     "shell.execute_reply.started": "2021-12-10T12:55:07.465483Z"
    },
    "papermill": {
     "duration": 0.024602,
     "end_time": "2022-01-02T18:59:44.412066",
     "exception": false,
     "start_time": "2022-01-02T18:59:44.387464",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def search_data(sequence_length, num_of_depend, label_start_idx,num_for_predict, units, points_per_hour):\n",
    "    '''\n",
    "    Parameters\n",
    "    ----------\n",
    "    sequence_length: int, length of all history data\n",
    "    num_of_depend: int,\n",
    "    label_start_idx: int, the first index of predicting target\n",
    "    num_for_predict: int, the number of points will be predicted for each sample\n",
    "    units: int, week: 7 * 24, day: 24, recent(hour): 1\n",
    "    points_per_hour: int, number of points per hour, depends on data\n",
    "    Returns\n",
    "    ----------\n",
    "    list[(start_idx, end_idx)]\n",
    "    '''\n",
    "\n",
    "    if points_per_hour < 0:\n",
    "        raise ValueError(\"points_per_hour should be greater than 0!\")\n",
    "\n",
    "    if label_start_idx + num_for_predict > sequence_length:\n",
    "        return None\n",
    "\n",
    "    x_idx = []\n",
    "    for i in range(1, num_of_depend + 1):\n",
    "        start_idx = label_start_idx - points_per_hour * units * i\n",
    "        end_idx = start_idx + num_for_predict\n",
    "        if start_idx >= 0:\n",
    "            x_idx.append((start_idx, end_idx))\n",
    "        else:\n",
    "            return None\n",
    "\n",
    "    if len(x_idx) != num_of_depend:\n",
    "        return None\n",
    "\n",
    "    return x_idx[::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:44.446912Z",
     "iopub.status.busy": "2022-01-02T18:59:44.446205Z",
     "iopub.status.idle": "2022-01-02T18:59:44.453515Z",
     "shell.execute_reply": "2022-01-02T18:59:44.453962Z"
    },
    "papermill": {
     "duration": 0.030037,
     "end_time": "2022-01-02T18:59:44.454113",
     "exception": false,
     "start_time": "2022-01-02T18:59:44.424076",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_sample_indices(data_sequence, num_of_weeks, num_of_days, num_of_hours, label_start_idx, num_for_predict, points_per_hour=12):\n",
    "    '''\n",
    "    Parameters\n",
    "    ----------\n",
    "    data_sequence: np.ndarray shape is (sequence_length, num_of_vertices, num_of_features)\n",
    "    num_of_weeks, num_of_days, num_of_hours: int\n",
    "    label_start_idx: int, the first index of predicting target\n",
    "    num_for_predict: int,the number of points will be predicted for each sample\n",
    "    points_per_hour: int, default 12, number of points per hour\n",
    "    Returns\n",
    "    ----------\n",
    "    week_sample: np.ndarray shape is (num_of_weeks * points_per_hour, num_of_vertices, num_of_features)\n",
    "    day_sample: np.ndarray shape is (num_of_days * points_per_hour,  num_of_vertices, num_of_features)\n",
    "    hour_sample: np.ndarray   shape is (num_of_hours * points_per_hour, num_of_vertices, num_of_features)\n",
    "    target: np.ndarray shape is (num_for_predict, num_of_vertices, num_of_features)\n",
    "    '''\n",
    "    week_sample, day_sample, hour_sample = None, None, None\n",
    "#------------------------------------Ignore\n",
    "    if label_start_idx + num_for_predict > data_sequence.shape[0]: \n",
    "        return week_sample, day_sample, hour_sample, None\n",
    "\n",
    "    if num_of_weeks > 0:\n",
    "        week_indices = search_data(data_sequence.shape[0], num_of_weeks, label_start_idx, num_for_predict,7 * 24, points_per_hour)\n",
    "        if not week_indices:\n",
    "            return None, None, None, None\n",
    "\n",
    "        week_sample = np.concatenate([data_sequence[i: j] for i, j in week_indices], axis=0)\n",
    "\n",
    "    if num_of_days > 0:\n",
    "        day_indices = search_data(data_sequence.shape[0], num_of_days,  label_start_idx, num_for_predict, 24, points_per_hour)\n",
    "        if not day_indices:\n",
    "            return None, None, None, None\n",
    "\n",
    "        day_sample = np.concatenate([data_sequence[i: j] for i, j in day_indices], axis=0)\n",
    "#----------------------------------Continue\n",
    "    if num_of_hours > 0:\n",
    "        hour_indices = search_data(data_sequence.shape[0], num_of_hours, label_start_idx, num_for_predict, 1, points_per_hour)\n",
    "        if not hour_indices:\n",
    "            return None, None, None, None\n",
    "        hour_sample = np.concatenate([data_sequence[i: j] for i, j in hour_indices], axis=0)\n",
    "    \n",
    "    if num_of_hours > 10:\n",
    "        return 1;\n",
    "    target = data_sequence[label_start_idx: label_start_idx + num_for_predict]\n",
    "\n",
    "    return week_sample, day_sample, hour_sample, target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:44.482208Z",
     "iopub.status.busy": "2022-01-02T18:59:44.481452Z",
     "iopub.status.idle": "2022-01-02T18:59:44.499450Z",
     "shell.execute_reply": "2022-01-02T18:59:44.499922Z",
     "shell.execute_reply.started": "2021-12-10T12:39:46.078756Z"
    },
    "papermill": {
     "duration": 0.033739,
     "end_time": "2022-01-02T18:59:44.500082",
     "exception": false,
     "start_time": "2022-01-02T18:59:44.466343",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def read_and_generate_dataset(graph_signal_matrix_filename, num_of_weeks, num_of_days, num_of_hours, num_for_predict, points_per_hour=12):\n",
    "    '''\n",
    "    Parameters\n",
    "    ----------\n",
    "    graph_signal_matrix_filename: str, path of graph signal matrix file\n",
    "    num_of_weeks, num_of_days, num_of_hours: int\n",
    "    num_for_predict: int\n",
    "    points_per_hour: int, default 12, depends on data\n",
    "    Returns\n",
    "    ----------\n",
    "    feature: np.ndarray, shape is (num_of_samples, num_of_depend * points_per_hour, num_of_vertices, num_of_features)\n",
    "    target: np.ndarray, shape is (num_of_samples, num_of_vertices, num_for_predict)\n",
    "    '''\n",
    "    #--------------------------------- Read original data \n",
    "    data_seq = np.load(graph_signal_matrix_filename)['data']  # (sequence_length, num_of_vertices, num_of_features) (16992, 307, 3)\n",
    "    \n",
    "    #---------------------------------\n",
    "    all_samples = []\n",
    "    for idx in range(data_seq.shape[0]):\n",
    "        sample = get_sample_indices(data_seq, num_of_weeks, num_of_days, num_of_hours, idx, num_for_predict, points_per_hour)\n",
    "        if ((sample[0] is None) and (sample[1] is None) and (sample[2] is None)):\n",
    "            continue\n",
    "\n",
    "        week_sample, day_sample, hour_sample, target = sample #  week_sample, day_sample are None because we are predicting per hour\n",
    "        #print(target.shape) # hour_sample and target (12, 307, 3)\n",
    "        sample = []  # [(week_sample),(day_sample),(hour_sample),target,time_sample]\n",
    "#-------------------------------- Ignore\n",
    "        if num_of_weeks > 0:\n",
    "            week_sample = np.expand_dims(week_sample, axis=0).transpose((0, 2, 3, 1))  # (1,N,F,T)\n",
    "            sample.append(week_sample)\n",
    "\n",
    "        if num_of_days > 0:\n",
    "            day_sample = np.expand_dims(day_sample, axis=0).transpose((0, 2, 3, 1))  # (1,N,F,T)\n",
    "            sample.append(day_sample)\n",
    "#----------------------------------Continue\n",
    "        if num_of_hours > 0:\n",
    "            hour_sample = np.expand_dims(hour_sample, axis=0).transpose((0, 2, 3, 1))  # (1,N,F,T)\n",
    "            sample.append(hour_sample)\n",
    "\n",
    "        target = np.expand_dims(target, axis=0).transpose((0, 2, 3, 1))[:, :, 0, :]  # (1,N,T)\n",
    "        sample.append(target)\n",
    "        time_sample = np.expand_dims(np.array([idx]), axis=0)  # (1,1)\n",
    "        sample.append(time_sample)\n",
    "        all_samples.append(sample)#sampe：[(week_sample),(day_sample),(hour_sample),target,time_sample] = [(1,N,F,Tw),(1,N,F,Td),(1,N,F,Th),(1,N,Tpre),(1,1)]\n",
    "\n",
    "    split_line1 = int(len(all_samples) * 0.6)\n",
    "    split_line2 = int(len(all_samples) * 0.8)\n",
    "\n",
    "    training_set = [np.concatenate(i, axis=0)  for i in zip(*all_samples[:split_line1])] #[(B,N,F,Tw),(B,N,F,Td),(B,N,F,Th),(B,N,Tpre),(B,1)]\n",
    "    validation_set = [np.concatenate(i, axis=0) for i in zip(*all_samples[split_line1: split_line2])]\n",
    "    testing_set = [np.concatenate(i, axis=0) for i in zip(*all_samples[split_line2:])]\n",
    "\n",
    "    return training_set, validation_set, testing_set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:44.527565Z",
     "iopub.status.busy": "2022-01-02T18:59:44.526797Z",
     "iopub.status.idle": "2022-01-02T18:59:46.159491Z",
     "shell.execute_reply": "2022-01-02T18:59:46.160177Z",
     "shell.execute_reply.started": "2021-12-10T12:39:46.98975Z"
    },
    "papermill": {
     "duration": 1.648016,
     "end_time": "2022-01-02T18:59:46.160324",
     "exception": false,
     "start_time": "2022-01-02T18:59:44.512308",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(16992, 307, 3)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "graph_signal_matrix_filename = '../input/pems-dataset/data/PEMS04/PEMS04.npz'\n",
    "data = np.load(graph_signal_matrix_filename)\n",
    "data['data'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:46.189644Z",
     "iopub.status.busy": "2022-01-02T18:59:46.188943Z",
     "iopub.status.idle": "2022-01-02T18:59:50.241384Z",
     "shell.execute_reply": "2022-01-02T18:59:50.241952Z",
     "shell.execute_reply.started": "2021-12-10T12:39:47.766013Z"
    },
    "papermill": {
     "duration": 4.068513,
     "end_time": "2022-01-02T18:59:50.242110",
     "exception": false,
     "start_time": "2022-01-02T18:59:46.173597",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "data = np.load(graph_signal_matrix_filename)\n",
    "data['data'].shape\n",
    "\n",
    "num_of_vertices = 307\n",
    "points_per_hour = 12\n",
    "num_for_predict = 12\n",
    "num_of_weeks = 0\n",
    "num_of_days = 0\n",
    "num_of_hours = 1\n",
    "\n",
    "training_set, validation_set, testing_set = read_and_generate_dataset(graph_signal_matrix_filename, 0, 0, num_of_hours, \n",
    "                                                                      num_for_predict, points_per_hour=points_per_hour)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.013082,
     "end_time": "2022-01-02T18:59:50.268271",
     "exception": false,
     "start_time": "2022-01-02T18:59:50.255189",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Data normalization\n",
    "\n",
    " the data are transformed by zero-mean normalization x′ = x −mean(x) to let the average be 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:50.299157Z",
     "iopub.status.busy": "2022-01-02T18:59:50.298455Z",
     "iopub.status.idle": "2022-01-02T18:59:50.307298Z",
     "shell.execute_reply": "2022-01-02T18:59:50.307778Z",
     "shell.execute_reply.started": "2021-12-10T13:03:05.96189Z"
    },
    "papermill": {
     "duration": 0.025887,
     "end_time": "2022-01-02T18:59:50.307951",
     "exception": false,
     "start_time": "2022-01-02T18:59:50.282064",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def normalization(train, val, test):\n",
    "    '''\n",
    "    Parameters\n",
    "    ----------\n",
    "    train, val, test: np.ndarray (B,N,F,T)\n",
    "    Returns\n",
    "    ----------\n",
    "    stats: dict, two keys: mean and std\n",
    "    train_norm, val_norm, test_norm: np.ndarray,\n",
    "                                     shape is the same as original\n",
    "    '''\n",
    "\n",
    "    assert train.shape[1:] == val.shape[1:] and val.shape[1:] == test.shape[1:]  # ensure the num of nodes is the same\n",
    "    mean = train.mean(axis=(0,1,3), keepdims=True)\n",
    "    std = train.std(axis=(0,1,3), keepdims=True)\n",
    "    print('mean.shape:',mean.shape)\n",
    "    print('std.shape:',std.shape)\n",
    "\n",
    "    def normalize(x):\n",
    "        return (x - mean) / std\n",
    "\n",
    "    train_norm = normalize(train)\n",
    "    val_norm = normalize(val)\n",
    "    test_norm = normalize(test)\n",
    "\n",
    "    return {'_mean': mean, '_std': std}, train_norm, val_norm, test_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:50.337613Z",
     "iopub.status.busy": "2022-01-02T18:59:50.336930Z",
     "iopub.status.idle": "2022-01-02T18:59:55.121356Z",
     "shell.execute_reply": "2022-01-02T18:59:55.120773Z"
    },
    "papermill": {
     "duration": 4.800485,
     "end_time": "2022-01-02T18:59:55.121487",
     "exception": false,
     "start_time": "2022-01-02T18:59:50.321002",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean.shape: (1, 1, 3, 1)\n",
      "std.shape: (1, 1, 3, 1)\n"
     ]
    }
   ],
   "source": [
    "train_x = np.concatenate(training_set[:-2], axis=-1)  # (B,N,F,T')\n",
    "val_x = np.concatenate(validation_set[:-2], axis=-1)\n",
    "test_x = np.concatenate(testing_set[:-2], axis=-1)\n",
    "\n",
    "train_target = training_set[-2]  # (B,N,T)\n",
    "val_target = validation_set[-2]\n",
    "test_target = testing_set[-2]\n",
    "\n",
    "train_timestamp = training_set[-1]  # (B,1)\n",
    "val_timestamp = validation_set[-1]\n",
    "test_timestamp = testing_set[-1]\n",
    "\n",
    "(stats, train_x_norm, val_x_norm, test_x_norm) = normalization(train_x, val_x, test_x)\n",
    "\n",
    "all_data = {'train': { 'x': train_x_norm, 'target': train_target,'timestamp': train_timestamp},\n",
    "            'val': {'x': val_x_norm, 'target': val_target, 'timestamp': val_timestamp},\n",
    "            'test': {'x': test_x_norm, 'target': test_target, 'timestamp': test_timestamp},\n",
    "            'stats': {'_mean': stats['_mean'], '_std': stats['_std']} }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:55.158986Z",
     "iopub.status.busy": "2022-01-02T18:59:55.158295Z",
     "iopub.status.idle": "2022-01-02T18:59:55.165069Z",
     "shell.execute_reply": "2022-01-02T18:59:55.166150Z",
     "shell.execute_reply.started": "2021-12-10T12:05:55.889468Z"
    },
    "papermill": {
     "duration": 0.031431,
     "end_time": "2022-01-02T18:59:55.166358",
     "exception": false,
     "start_time": "2022-01-02T18:59:55.134927",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train x: (10181, 307, 3, 12)\n",
      "train target: (10181, 307, 12)\n",
      "train timestamp: (10181, 1)\n",
      "\n",
      "val x: (3394, 307, 3, 12)\n",
      "val target: (3394, 307, 12)\n",
      "val timestamp: (3394, 1)\n",
      "\n",
      "test x: (3394, 307, 3, 12)\n",
      "test target: (3394, 307, 12)\n",
      "test timestamp: (3394, 1)\n",
      "\n",
      "train data _mean : (1, 1, 3, 1) [[[[2.07227338e+02]\n",
      "   [5.13195612e-02]\n",
      "   [6.34740574e+01]]]]\n",
      "train data _std : (1, 1, 3, 1) [[[[1.56477655e+02]\n",
      "   [4.78541626e-02]\n",
      "   [8.10351724e+00]]]]\n"
     ]
    }
   ],
   "source": [
    "print('train x:', all_data['train']['x'].shape)\n",
    "print('train target:', all_data['train']['target'].shape)\n",
    "print('train timestamp:', all_data['train']['timestamp'].shape)\n",
    "print()\n",
    "print('val x:', all_data['val']['x'].shape)\n",
    "print('val target:', all_data['val']['target'].shape)\n",
    "print('val timestamp:', all_data['val']['timestamp'].shape)\n",
    "print()\n",
    "print('test x:', all_data['test']['x'].shape)\n",
    "print('test target:', all_data['test']['target'].shape)\n",
    "print('test timestamp:', all_data['test']['timestamp'].shape)\n",
    "print()\n",
    "print('train data _mean :', all_data['stats']['_mean'].shape, all_data['stats']['_mean'])\n",
    "print('train data _std :', all_data['stats']['_std'].shape, all_data['stats']['_std'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-01-02T18:59:55.204554Z",
     "iopub.status.busy": "2022-01-02T18:59:55.202501Z",
     "iopub.status.idle": "2022-01-02T19:00:42.325357Z",
     "shell.execute_reply": "2022-01-02T19:00:42.326407Z",
     "shell.execute_reply.started": "2021-12-10T12:05:57.459518Z"
    },
    "papermill": {
     "duration": 47.145799,
     "end_time": "2022-01-02T19:00:42.326581",
     "exception": false,
     "start_time": "2022-01-02T18:59:55.180782",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "save file: ./PEMS04_r1_d0_w0_astcgn\n"
     ]
    }
   ],
   "source": [
    "file = os.path.basename(graph_signal_matrix_filename).split('.')[0]\n",
    "dirpath = '.'\n",
    "filename = os.path.join(dirpath, file + '_r' + str(num_of_hours) + '_d' + str(num_of_days) + '_w' + str(num_of_weeks)) + '_astcgn'\n",
    "print('save file:', filename)\n",
    "np.savez_compressed(filename,\n",
    "                train_x=all_data['train']['x'],train_target=all_data['train']['target'],train_timestamp=all_data['train']['timestamp'],\n",
    "                val_x=all_data['val']['x'], val_target=all_data['val']['target'],val_timestamp=all_data['val']['timestamp'],\n",
    "                test_x=all_data['test']['x'], test_target=all_data['test']['target'], test_timestamp=all_data['test']['timestamp'],\n",
    "                mean=all_data['stats']['_mean'], std=all_data['stats']['_std'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.549297,
     "end_time": "2022-01-02T19:00:42.890522",
     "exception": false,
     "start_time": "2022-01-02T19:00:42.341225",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "papermill": {
   "duration": 63.803621,
   "end_time": "2022-01-02T19:00:43.947524",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2022-01-02T18:59:40.143903",
   "version": "2.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
